Universality classes and crossover behaviors in non-Abelian directed sandpiles 
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We study universality classes and crossover behaviors in non-Abelian directed sandpile models, 
in terms of the metastable pattern analysis. The non-Abelian property induces spatially correlated 
metastable patterns, characterized by the algebraic decay of the grain density along the propagation 
direction of an avalanche. Crossover scaling behaviors are observed in the grain density due to the 
interplay between the toppling randomness and the parity of the threshold value. In the presence of 
such crossovers, we show that the broadness of the grain distribution plays a crucial role in resolving 
the ambiguity of the universality class. Finally, we claim that the metastable pattern analysis is 
important as much as the conventional analysis of avalanche dynamics. 

PACS numbers: 05.70.Ln, 05.65.+b, 64.60.Ht 
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I. INTRODUCTION 

Scale invariance in avalanche systems is ubiquitously 
observed in nature, such as earthquakes, sandpiles, and 
Barkhausen avalanches of ferromagnetic materials. The 
statistics of avalanches follows some universal power- 
law distribution Metastable patterns between 
avalanches also exhibit spatially long-range correlations. 
A fractal structure in the crust of the earth formed by 
seismic events is one of such examples. It is, however, still 
questionable to find generic mechanisms for avalanche 
systems that can explain both the statistics of avalanches 
and long-range correlations in metastable patterns. 

To reveal the underlying common mechanism of scale 
invariance in avalanche systems, Bak, Tang, and Wiesen- 
feld (BTW) first proposed the paradigm of self-organized 
criticality with an Abelian deterministic version of undi- 
rected sandpile models. Since then, lots of BTW variants 
have been suggested and studied [4- 6] . In modeling sand- 
piles, the balance between slow driving and dissipation 
is the key ingredient. Grains are slowly added, toppled 
instantly whenever the instability threshold is overcome, 
and finally dissipated at boundaries of the system. It has 
been tested whether the universality class of avalanche 
dynamics can be changed by the breaking of Abelian 
symmetry [tJ or the consideration of stochasticity [8[ in 
toppling rules. The issue is still controversial due to con- 
flicting numerical results [1-G3]. The Abelian symmetry 
here means that the order of toppling events does not 
affect the final state. 

Contrary to undirected sandpile models, directed 
sandpile models (DSMs) with a preferred direction of 
avalanche propagation turn out to be more tractable ana- 
lytically as long as they have the Abelian symmetry [T^ — 
KL7| . This is because the Abelian symmetry in DSMs lets 
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metastable grain patterns be fully uncorrelated. How- 
ever, most real avalanche systems often exhibit spatially 
correlated grain patterns. It means that Abelian DSMs 
are not suitable enough to describe such systems and 
non-Abelian DSMs are more natural to be considered. 

When the Abelian symmetry is broken in DSMs [Tcl - 
[20j . spatially long-range correlations emerge in 
metastable grain and scar patterns. Scars are de- 
fined as traces of avalanche boundaries. Both scar and 
grain densities decay algebraically along the preferred 
direction of avalanche propagation with the same decay 
exponent. In our earlier work [20j |. additional scaling 
relations arc derived from avalanche flow equations. 
All avalanche exponents can be written in terms of the 
grain density exponent. With the stochastic toppling 
rule, the broken Abelian symmetry seems not to change 
the universality class. However, it changes the scaling 
behavior of metastable patterns with the non-zero value 
of the grain density exponent. When the toppling rule 
becomes deterministic, the Abelian symmetry turns out 
to be relevant to the universality class and non-Abelian 
DSMs belong to the mean-field (MF) universality class. 
These results were numerically confirmed in the simplest 
models [13]. 

However, the simplest discrete DSM with the non- 
Abelian deterministic toppling rule cannot be unique due 
to the interplay between the toppling rule and the lattice 
structure. For an example, if the number of grains at an 
unstable site is 3 and the number of its neighboring sites 
is 2, there are several ways to topple the last grain into 
one of the neighboring sites. Two versions of the top- 
pling bias for the last grain were considered, which lead 
to some deviation from MF results [20| . 

In this paper, we check out the validity of universal- 
ity classes in non-Abelian DSMs against the change of 
the toppling bias for the last grain and the instability 
threshold value. We also report various crossover behav- 
iors between universality classes. One might expect that 
the toppling bias effect becomes negligible as the thresh- 
old value increases. It turns out to be only true for the 
fully biased case. For alternatively and/or partially bi- 



ased cases, the parity of the threshold value plays a cru- 
cial role in avalanche dynamics. If the threshold value is 
odd, the stochasticity in toppling rules becomes relevant. 
As a result, the deterministic models belong to the non- 
Abelian stochastic universality class. It is somehow sur- 
prising because one believes the threshold value in sand- 
pile models to be irrelevant to universality classes, except 
for crossover behaviors discussed in undirected cases by 
Liibeck HU. 

This paper is organized as follows. In SecHIl we br iefly 
describe various DSMs including those proposed in [20j . 
In Sec. IIII1 we discuss the role of the threshold value 
and the toppling bias in determining universality classes 
by the conventional analysis of avalanche dynamics and 
the metastable pattern analysis. We also argue possible 
scenarios for universality classes and crossover behaviors, 
which are numerically confirmed in Sec. IIV1 Finally, we 
conclude this paper in Sec. [V] with the summary of main 
results and some suggestions. 



II. DIRECTED SANDPILE MODELS 



We consider DSMs defined on a two-dimensional \- 
rotated lattice of sizes (L±, L»). For convenience, we call 
L± and L|| as L and T, respectively. The preferred direc- 
tion of avalanche propagation is denoted by the "layer" 
t = 0,1, ...,T — 1 with open boundary conditions. The 
transverse direction is denoted by i — 0,1,--- ,L — 1 
with periodic boundary conditions (see Fig. Q}. Initially, 
to each site of the lattice an integer value (the number 
of grains), z(i,t) g {0, 1, z c — 1}, is assigned. Here z c 
denotes the threshold value. 

Given a stable configuration where z(i,t) < z c for all 
sites (i, t), one grain is added at a randomly chosen site on 
the top layer, z(i, 0) — > z(i, 0) + 1. For any unstable site 
with z{i,t) > z c , including that on the top layer, grains 
topple down to its left and right nearest-neighboring sites 



(a) 




FIG. 



t=T-l 



1: A two-dimensional ^-rotated lattice for DSMs. 




FIG. 2: Toppling rules of non-Abelian deterministic DSMs 
on a two-dimensional tilted lattice 20]: (a) aND, (b) bND, 
and (c) cND. The gray colored grains and arrow represent the 
state before toppling and the black colored ones represent the 
state after toppling. 



on the next layer t + 1 . 
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Here Ai jt represents the number of grains outgoing from 
(i,t) and Ai±i it+ i does the number of grains incoming 
to (i ± 1, t+ 1) from (i, t). The number of toppled grains 
is locally conserved: A M = A;_i it+1 + A. i+M+:L . 

The toppling at one layer may cause another toppling 
on the next layer. At unstable sites on the bottom layer 
t = T—l, toppled grains are dissipated out of the system. 
Only after another stable configuration is recovered by a 
series of toppling events, denoting an avalanche, a new 
grain is added on the top layer, t = 0, to keep generating 
another avalanche. 

By setting the toppling matrices, A and A, one may 
consider several variants of DSMs. Without losing any 
generality, we set A^ jt = z c in the Abelian version and 
A^.t = z(i,t) > z c in the non-Abelian version. The val- 
ues of Ai±i jt+ i can be either stochastically or determin- 
istically determined. For the stochastic case, each grain 
at the unstable site is toppled at random to one of the 
nearest neighbors at the next layer. Abelian and non- 
Abelian stochastic DSMs are denoted as AS and NS for 
short, respectively. For Abelian deterministic DSMs (AD 
for short), one should fix the values of Ai±i jt+ i under the 
condition Ai_i it _|_i + A.; + i it+ i = z c . For an example, each 
value of Ai-ti it+ i is set as one half of z c when z c is even. 

For non-Abelian deterministic DSMs (ND for short) 
where the simplest setup is not clear, the following three 
versions are considered in 1201: 



(i) A 



j±i,t+i 



(ii) A i±M+ i 

(iii) A i±M+ i 



k iiz(i,t)=2k 
k + £j±i,a(i,t) if z(i,t) = 2k + 1, 

k if z(i,t) = 2k 

k + <5i±i,i+i if z(i, t) = 2k + 1, 

z(i,t)/2. (2) 



Here, k is a positive integer and Sij denotes the Kronecker 
delta function. We call (i) the alternative bias version 
(aND) [l9[, (ii) the full bias version (bND), and (iii) the 
continuous version without bias (cND). For the aND, a 




"toppling arrow" a(i,t) of each site initially points to one 
of its neighbors, say, i — 1, in Fig. [2ja) . Whenever each 
grain is toppled at the unstable site, the direction of the 
arrow flips to the other neighbor. Figure [2] shows the case 
of z(i,t) = 3 and z c = 2. It is worthwhile to note that 
any toppling bias here is applied only to the last grain 
when the number of toppled grains at the unstable site, 
z(i,t), is odd. 

In order to check the validity of earlier ND versions, we 
introduce another ND version that covers both the aND 
and the bND by controlling the bias of the last grain at 
the toppled site as follows: 

A _(k if z(i,t) = 2k 

(iv) A J±M+1 -| fc + 5i±lr{it) if z[h t ) = 2k+ 1. 

(3) 

Here r(i, t) takes i + 1 with a probability p or i — 1 with a 
probability 1— p for each toppling event. For convenience, 
only the range of \ < p < 1 is considered. We call (iv) 
the partial bias version (pND for short). This model is 
exactly the same as the bND when p = 1. Moreover, 
it is almost the same as the aND when p = i, except 
that toppling arrows are initially quenched in the aND 
or completely annealed in the pND. 

An avalanche of DSMs can be quantified as mass s (the 
total number of toppled grains), duration t (the number 
of affected layers), area a (the total number of distinct 
toppled sites), width w (the distance between avalanche 
boundaries), and height h (the number of toppled grains 
per toppled site). There is no characteristic scale in 
avalanche dynamics, except for T as long as L is suffi- 
ciently larger than the maximum width of the avalanche. 
It is well known that all probability distribution func- 
tions of avalanche quantities are expected to follow the 



simple scaling form as 

P(x,T) ~ x- T *f (x/T D *) (4) 

for x £ {s, t, a, w, h}. Hence, any two quantities x and y 
scale, using the conditional expectation value, as 

E\y\x] ~ x 1 ^. (5) 

The relative exponent "f yx = (t x — 1)/(t v — 1) = D y /D x 
is derived from the identity P(x 7 T)dx — P{y,T)dy. 

We take full advantage of the relations, D t — 1 and 
(s) ~ T in DSMs, where (•) denotes the ensemble average. 
Based on the reasonable assumption for the compactness 
of the avalanche in low dimensional systems, i.e., a ^ wt 
and s ~ ah, we obtain the following scaling relations: 

Ixt = Dx for any x, 

D,(2-t.) = 1, (6) 
D a = D w + 1, and D s = D a + D h , 

which implies only two independent exponents left. 

Concerning metastable patterns of DSMs, we are able 
to define two additional scaling exponents. In order to 
quantify metastable patterns, we denote the probability 
to find z remaining grains at each layer t as 

Q(z,t)=Pr[z(i,t) = z}, (7) 

where Y^z^o Q( z > *) = 1- Along the avalanche propaga- 
tion direction, one can measure the grain density as 

L-l z c -l 

p(i) = z 5>(M)}= $>Q(M)' (8) 

i=0 z=0 

where (•) denotes the recurrent configuration average. 
Similarly, the fraction of occupied sites at each layer, 
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namely, the grain occupation density, can be measured 
as 



L-l 



P s (t) 



^X>[*(m)]> = Eg(m)' 



(9) 



i=0 



where 8(z) denotes the Heaviside step function that gives 
for z = and 1 for z > 0. Both p(t) and p g (t) decay 
algebraically with the grain density exponent a at large 
t in non-Abelian DSMs as 



(10) 



For Abclian DSMs with an uniform grain density (a = 0), 
we consider the scar density p sc (t) and the scar exponent 
a sc with the same definition of pit) where z(i,t) is just 
replaced by the quantity, b(i,t). The value of b(i,t) is 
1 if the site (i,t) has recently been a part of avalanche 
boundaries, or otherwise. We call {b(i,t)} a scar. 



III. AVALANCHE DYNAMICS 

We begin with the discussion of statistical properties 
observed in metastable grain and scar patterns. The 
metastable pattern analysis is important as much as 
the conventional analysis of avalanche dynamics. As re- 
ported in [20l | . the scar exponent is a good indicator for 
the universality classes of DSMs since it is directly related 
to the exponent of the avalanche width as a sc = D w . It 
is because the scar density is inversely proportional to 
the typical avalanche width. 



L-l 



p sc {t) 



£5>(M)) 



1 L 



i=0 



Lw(t) 



w(t)-\ 



(11) 



It is worthwhile to note that grains in non-Abelian DSMs 
can remain only at avalanche boundaries, which implies 
a = a sc . On the other hand, Abelian DSMs exhibit 
an uncorrelated and uniform grain density in metastable 
grain patterns, i.e., a — 0, irrespective of the a sc value. 

We intuitively argue the interplay between avalanche 
flow and metastable grain or scar patterns in DSMs. Let 
us define N(t) as the number of grains toppled from the 
layer t to the next layer t+1 within an avalanche. It scales 
as N(t) ~ w(t)h(t) ~ t D w+Dh w ith the assumption for 
the compactness of an avalanche. The evolution of N(t), 
avalanche flow, can be written as 



dN(t) 
dt 



N(t) - N(t- 1) 



iew(t) 



(12) 



where n(i,t) denotes the amount of avalanche flow at the 
site (i,t). The value of n(i,t) depends on both the value 
of z(i, t) before toppling and the number of grains toppled 
from the previous layer to the site (i, t). The summation 
is over avalanche boundaries and bulk sites, belonging to 
w(t). Here avalanche boundaries denote the outermost 
sites of an avalanche at a given layer and avalanche bulk 
sites do the sites between avalanche boundaries. 



A. Abelian cases 

In AD DSMs, it is well known that metastable grain 
patterns are fully uncorrelated and that probabilities to 
find z(i,t) € {0, 1, z c — 1} grains at a site are the same 
as 1/ z c |l4[ . The number of toppled grains at every un- 
stable site is z c , such that D h = 0. Since n(i,t) = 
at avalanche bulk sites, only avalanche boundaries con- 
tribute to avalanche flow. At the right avalanche bound- 
ary, the avalanche width w increases by one if the top- 
pling at the rightmost site of the avalanche, say (i,t), 
leads to a toppling at a site (i + l,t + l). The probability 
is given by 



A i+ i 



t+i 



Pr[*(i + l,t+l) + A i+1 , t+1 > z c ] 



The probability to decrease w by 1 is 1 IP ah ■ Likewise, 
probabilities to increase and decrease w by one at the left 
avalanche boundary are 1 — p AD and p Ar> , respectively. 
For the range of < p AB < 1, avalanche boundaries 
in the AD are exactly mapped onto the spatiotemporal 
trajectories of random walks starting at the same site. 
The avalanche flow equation is therefore described as 



dN(t) 
dt 



v(t), 



(13) 



where T)(t) denotes an uncorrelated noise with zero mean 
and unit variance. As a result, we obtain a sc = D w = | 
with Dh — 0, which yields all avalanche exponent values. 
We note that the result in the AD is robust, independent 
of the value of z c and any possible toppling bias in the 
case of p AD ^ §. 

In AS DSMs, the fluctuations of Aj±i.t+i result- 
ing from the stochastic toppling rule enable multiple 
topplings (more than one toppling at the same site). 
Metastable grain patterns are still fully uncorrelated 
because of the characteristic of the Abelian symmetry 
(a = 0) [H, In that sense, avalanche boundaries in 
the AS are also mapped onto the trajectories of random 
walks, a sc = D w — |. However, avalanche flow in the 
AS is mainly affected by avalanche bulk sites. In gen- 
eral, the numbers of grains that bulk sites receive from 
the previous layer can be written as mz c + 1, where m is a 
non-negative integer and < I < z c — 1. If z(i, t) + l < z c 
with the value of z(i,t) in the metastable state, / grains 
would be left behind avalanche flow, i.e., n(i,t) = —I. If 
z(i, t) + 1 > z c , avalanche flow would sweep z c — I grains, 
i.e., n(i, t) = z c — I. Based on the reasonable assumption 
that the probability distribution of I is uniform, n(i,t) is 
considered as a random variable r)(i,t) with zero mean 
and finite variance. As a result, the avalanche flow equa- 
tion in the AS is described as 



dN(t) 
dt 



iew(t) 



VW)v(t), (14) 



which gives Dh — j- A positive value of Dh reflects the 
existence of the multiple toppling. Once again, above 



5 



results are also robust, independent of the value of z c . 
The toppling bias can be considered as setting the prob- 
abilities of each grain toppled to the right and the left 
nearest neighbors at the next layer with p AS and 1 — p AS , 
respectively. In other words, 



Pa 



(A 



i+\,t+l) 



where (•) denotes the ensemble average due to the 
stochasticity in toppling rules. When < p AS < 1, 
such toppling bias does not change the scaling proper- 
ties of avalanche dynamics (numerically confirmed but 
not shown here). 



B. Non-Abelian cases 

In NS DSMs, spatially correlated metastable patterns 
are well described as the algebraic decay of the grai n den- 
sity along the avalanche propagation direction [18| . Due 
to such spatial correlations in metastable patterns, the 
uncorrelated noise rj in the avalanche flow equation for 
Abelian cases should be replaced by some correlated noise 
to represent NS metastable patterns. Naively speaking, 
one can interpret 77 in Eq. (|14p as p sc having the same 
dimension of i -1 / 2 as rj. Thus, we replace p sc by p in 
the NS, while keeping \/w(t) due to the stochasticity in 
toppling rules: 



(15) 



which yields Dh — 1 — |a. As a result all other avalanche 
exponents in the NS can be written in terms of a, such 
as t s = 2( 4 3 T° ) and r t = 2 - § . 

The value of a in the NS is numerically observed as 
0.45, slightly less than |. In the earlier work [20| . it 
is argued that the deviation seems to be attributed to 
logarithmic corrections to scaling. So we simply assume 
the random walk behaviors of avalanche boundaries, such 
that a = o) irrespective of the toppling bias and the 
threshold value (numerically confirmed but not shown 
here). Then, all avalanche exponents in the NS becomes 
exactly the same as those in the AS: r s = tj = |, 
and Dh = |. From now on, we call the set of these 
exponents as the "NS" universality class. One can say 
that the NS belongs to the same universality class of the 
AS in the following sense: for Abelian cases, avalanche 
flow can sweep and lose many grains at the same time 
due to the uniform grain density. On the other hand, 
for non- Abelian cases, avalanche flow can sweep only few 
grains due to the power-law decaying grain density, but it 
can also leave only few grains behind as taking all grains 
at unstable sites. For both cases, the scaling property 
of N(t) is apparently unaffected by the grain density 
as long as the toppling rule is stochastic. Finally, we 
naively propose another possible scenario for a = i in 



TABLE I: Avalanche exponents {r s , r t , Dh}, grain density 
and scar exponents, a and a sc and our conjecture in two- 
dimensional DSMs: the error bars of numerical results are 
obtained from effective exponent plots by the assumption of 
the dominant simple power-law scaling [23|]. It is worthwhile 
to note that D 3 = rt and Dt = 1 in all DSMs and a sc = a in 
non-Abelian DSMs. 



Model 



n 



D h 



Mean field (d u = 2) 

Abelian 

Deterministic 
Stochastic 

Non-Abelian 
Stochastic 

- z c = 2 

Deterministic 
aND 

- z c = 2 

- z c = 3 

- z c = 4 
bND 

- z c = 2 

- z c = 3 

- z c = 4 
cND 

- 2c = 2 



0, i 

0, § 



2(3-a) 
4— a 

1.43(1) 



a 



cy a 1 3a 

2 2 

1.78(1) 0.31(3) 0.45(3) 



1.49(1) 
1.46(2) 
1.53(1) 

1.43(1) 
1.46(2) 
1.51(2) 



1.94(1) 0.07(1) 0.86(3) 

1.78(3) 0.26(1) 0.37(1) 

2.00(5) 0.16(2) 0.97(1) 

1.79(1) 0.06(1) 0.69(5) 

1.85(3) 0.20(2) 0.78(14) 

2.02(3) 0.15(1) 0.89(15) 



1.52(3) 2.04(3) 0.06(4) 0.91(11) 



the NS by mapping metastable patterns onto the (1+1)- 
dimensional spatiotemporal configuration of 2A — > A 
coagulation-diffusion model, where the particle density 
temporally decays as t~ x l 2 \2^. However, this scenario 
requires further investigation. 

In ND DSMs, there are no fluctuations if we focus 
on the continuous version, namely the cND, shown in 
Fig. [2jc). The avalanche flow equation can be written as 



dN(t) 
dt 



,(t)p(t), 



(16) 



yielding Dh = 1 — a. We find that the scaling expo- 
nents of avalanche mass and duration are MF values, in- 
dependent of a, i.e., t s — I and r t — 2, whereas other 
exponents still depend on a. In that sense, we point 
out that one should check other avalanche exponents be- 
sides those of the avalanche mass and duration to dis- 
cuss the universality class in the ND. If the assumption 
of a = D w = 1 is accepted based on the linear behav- 
ior of avalanche boundaries, all avalanche exponents turn 
to be MF values, except for the case of the avalanche 
width [24j . Furthermore, this may also correspond to 
the MF behavior of the (d+ l)-dimcnsional spatiotempo- 
ral configuration of the coagulation-diffusion model when 
d > d u = 2, where the particle density decays as i _1 [22J. 
The peculiar "MF" behavior in the two-dimensional ND, 
i.e., appeared in the low dimensional system, can be eas- 
ily understood in the context of the shape of N(t) with 
its width and height. 
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Toppling rules we considered in the ND suppress the 
fluctuations of the height profile in avalanche flow much 
more than those in the NS do. This leads to spread 
grains wider and wider and makes avalanche boundaries 
grow faster, almost ballistically. Such a positive feedback 
enables D w = 1 to be larger than \ for all other DSMs. 
The resultant Dh = indicates the MF behavior with 
a = 1, whereas Dh = for any dimensions in the AD. 

Table U presents the avalanche exponents, the grain 
density exponent, and the scar exponent, as well as the 
exact results for the Abelian case and our conjecture for 
the non- Abelian case. 

We report that the numerical results in the aND and 
the bND clearly deviate from MF results. As the possible 
origin of anomalous scaling behaviors, we claim that they 
are attributed to the toppling bias applied to the last 
grain when z(t, t) — 2k + 1 > z c in Eq. ©. The toppling 
bias effect in the ND is numerically tested as z c increases. 
We find that the increment of z c may yield either that 
the toppling bias becomes negligible compared to other 
(2k) grains, or that it turns to be relevant to change the 
universality class. The aND with even z c values and the 
bND with all z c values correspond to the former. They 
show the MF behavior with some logarithmic corrections 
for small z c values and without such corrections for large 
z c values as expected (see Table QJ. Scaling behaviors in 
the aND with odd z c values are clearly deviated from MF 
results, and seem to be those of the NS class. 

Two different mechanisms can be illustrated for the 
MF class and the NS class emerged in the aND against 
the parity of z c . After the transient period of sandpilc 
models, each site having ever toppled should be empty or 
occupied by at least [z c /2] grains. Here [m] is the integer 
value that is smaller than or equal to m. 

In the even case with z c = 2k, at least k grains at 
the rightmost toppled site of the avalanche, say, (i,t), 
topple to its right neighboring site on the next layer, 
(i + 1, t + 1). The site (i + 1, t + 1) will topple if occu- 
pied by grains, or will not if empty. At small t, the grain 
density is relatively high and the avalanche width is still 
small. Grains swept by avalanche flow quickly diffuse 
to avalanche boundaries, which causes the ballistic be- 
havior. On the other hand, at large t, remaining grains 
are rare in metastable patterns. Avalanche flow keeps 
losing grains at avalanche boundaries, which also causes 
the ballistic behavior, resulting in a = 1. 

In the odd case with z c = 2k + 1, the minimal number 
of grains at a stable site is k unless empty. Assume that 
z c grains are at the rightmost site of an avalanche, say 
(i, t), and its right neighboring site (i + 1, t + 1) is occu- 
pied by k grains. The site (i + 1, t + 1) will topple only 
if it receives k + 1 grains from the site (i,t), or will not 
if it receives k grains. Whether the site topples or not 
depends on the value of the toppling arrow a(i,t). This 
randomness inherent to toppling arrows leads to the ran- 
dom walk behaviors of avalanche boundaries, i.e., a = ^. 
If more than z c grains topple at the rightmost site or 
if its right neighboring site is occupied by more than k 



grains, the randomness of the alternative toppling bias 
is suppressed as in the aND with even z c values. The 
same scaling behavior is also observed in the pND with 
p = i, where the toppling bias is completely random. 
In the bND with no intrinsic randomness (the pND with 
p = 1), we do not expect the random walk behaviors 
of avalanche boundaries, irrespective of the parity of the 
threshold value. The pND in the range of \ < P < 1 will 
be discussed with numerical results in Sec. llVl 



IV. NUMERICAL RESULTS 

We performed extensive numerical simulations for all 
DSMs to confirm our conjecture for avalanche exponents 
in terms of the scar exponent, a sc (a — a sc for non- 
Abelian cases), up to L = T = 2 13 (T = 2 15 or L = 2T 
in some cases). 

Using conventional successive slope techniques for at 
most 10 9 avalanches in the steady state enough after 
the transient period, we measure the avalanche exponent 
set {t x ,D x } for all x and the grain density exponent a 
(a sc for Abelian cases). As discussed, spatially correlated 
scars are observed in all DSMs with the non-zero values of 
a sc . The spatial correlations of metastable patterns are 
only observed in non-Abelian DSMs with the non-zero 
values of a (see Tableland Fig. [3]). 

We calculate effective exponents from power-law distri- 
butions of various avalanche quantities, scaling relations 
among those quantities, and the algebraic decay of the 
grain occupation density in terms of following definitions: 



Df(t) 
a cS (t) 



lnP(x) 


— \nP(x/m) 


lnx - 


- \w(x/m) 


\nE[x\t] - 


- ]nE[x\t/m] 


kit- 


ln(i/m) 


ing (t) 


- \np e (t/m) 


Int- 


- ln(i/m) 



(17) 
(18) 
(19) 



as taking appropriate values of m. In order to get the 
reasonable resolution of effective exponents, we set m = 5 
for t cS and a cS and m — 1 for D cS , respectively. The 
exponents, t x and D x (= j x t), are directly obtained by 
the linear fit of scaling regimes in power-law distribution 
functions and scaling relations. Scaling regimes can be 
identified cither from the flat region of effective exponents 
or the systematic tendency as the system size increases. 
Regarding j yx , instead of the original definition, i?[y|a;] ~ 
x lyx , we make use of the modified definition, [ z/ 1 a;] = 
x + Cx' ty!C , when y > x for each avalanche 25]. Due 
to s > z c t and h > z c by definition, we obtain 7 st and 
7fc f using E[s\t] = z c t + Ci^* and E[h\t] = z c + C'f ht , 
respectively. The values of scaling exponents are listed 
in Table U which confirm our conjecture by means of the 
avalanche flow analysis in Sec. IIIII 

In the aND, the scaling behaviors of grain occupation 
densities definitely depend on the parity of the thresh- 
old value, z c , as shown in Fig. Ufa). For the even z c 
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FIG. 4: (Color online) Numerical results of the aND for various z c values: (a) grain occupation densities p e (t) versus !ona 
lattice with T = 2 14 for z c = 3, 5, 7, and 9 (T = 2 13 otherwise) and L = 2T. (b) Finite-size scaling of P(t, T) with the MF 
value (jt = 2) for z c = 2 (lines) and z c = 4 (open symbols), and with the NS value (r t = 7/4) at z c = 3 (filled symbols). Here, 
T = L — 2 10 , 2 11 , 2 12 , and 2 13 from right to left and D t = 1 for all cases, (c) Crossover scaling of grain occupation densities 
for odd z c values by p e (t)Zc as a function of t/z%, which are the same data as (a), (d) Rescaled grain distributions Q r (z,t) at 
two different layers, t = 16 (open symbols) and 4096 (filled symbols) for z c — 50 (A) and 51 (o) on a lattice with T = 2 13 and 
L = 2T. 
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FIG. 5: (Color online) (a) Grain occupation densities p g (t) versus t in the pND at a fixed z c = 3 as p increases from 0.5 (the 
aND as the thin line at the top) to 1 (the bND as the thick line at the bottom). We used a lattice with T = 2 13 and L = 2T. 
(b) Crossover scaling of grain occupation densities at z c = 3 by p e (t) /(l — p) 1 ^ 2 as a function of t(l — p) 1//2 . 



values, the scaling regime of a cS (t) approaches 1 as z c 
increases (not shown) . The value of a at z c = 3 appears 
to be h within errors rather than 1. Interestingly, for 
the larger odd z c values, we observe the crossover be- 
havior from a cS (t) = 1 (with logarithmic corrections) to 
h at some "crossover layer," t x , which clearly increases 
as z c increases. We investigate the same parity effect 



on the universality class in terms of other effective ex- 
ponents, r° ff , T t cff , and D'jf , for avalanche distributions 
(not shown). The effective exponents for z c = 2 and 4 
approach MF values (r s = |, Tt = 2, and Dh = 0) or at 
least have the tendency to approach them. On the other 
hand, effective exponents at z c — 3 clearly deviate from 
MF values. We also confirm the finite-size scaling (FSS) 
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FIG. 6: (Color online) Numerical results of the bND for vari- 
ous z c values: (a) grain occupation densities p g (t) versus t on 
a lattice with T = 2 14 and L = 2T (T = 2 13 for z c = 50 and 
51), where lines for z c — 2, 4, and 50 are from top to bottom, 
(b) Finite-size scaling of P(t, T) with MF values for z c — 2 
(lines), z c — 3 (filled symbols), and z c — 4 (open symbols). 
Here T — L — 2 10 , 2 U , 2 12 , and 2 13 from right to left, (c) 
Rescaled grain distributions Q r (z,t) at two different layers, 
i = 16 (open symbols) and 4096 (filled symbols) for z c — 50 
(A) and 51 (o) on a lattice with T = 2 13 and L = 2T. 



behaviors of avalanche distributions with MF values for 
z c = 2 and 4, and with NS values (t s — r t — j, and 
Dh = j) at z c = 3. FigureSJb) shows the successful FSS 
collapses of avalanche duration distributions. 

Numerical data in the aND are averaged over a lot of 
avalanches (up to 10 9 with only few numbers of ensem- 
ble for random initial configuration setups) in the steady 
state after the transient. Since toppling arrows, {a(i,i)}, 
are annealed enough by the random addition of grains in 
the steady state, it is guaranteed that all the results are 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

z/(z c -l) 

FIG. 7: (Color online) (a) Grain occupation densities p g (t) 
versus t in the NS for various z c values (lines for z c — 2, 
20, 50, and 100 from top to bottom at the right). We used 
a lattice with T = 2 1S and L = 2T. (b) In particular to 
the cases for z c = 25 and 100, rescaled grain distributions, 
Q r {z,t), are measured at two different layers, t = 16 (open 
symbols) and 4096 (filled symbols). 



completely independent of initial conditions, which are 
also numerically confirmed. 

Figure @Jc) shows the crossover scaling of grain occu- 
pation densities in the aND for odd z c values. It can be 
written as 



p e (t)=t- 1 g(t/t x ), 



(20) 



where the scaling function behaves as g(x) ~ O(l) for 
small x and g(x) ~ \fx for large x. At small x (t < £ x ), 
grain occupation densities are relatively high while the 
avalanche width is still small. Grains at avalanche bulk 
sites quickly diffuse to avalanche boundaries. The num- 
bers of grains at avalanche boundaries tend to be larger 
than z c and thus make ballistic avalanche boundaries 
(a = 1). As the grain occupation densities decay to 
become low at large x (t > t x ), the numbers of grains at 
avalanche boundaries tend to be z c . Then the stochas- 
ticity of avalanche boundaries is disclosed (a = |). The 
value of t x indicates the crossover from the dense region 
to the sparse region of remaining grains. We present some 
naive test for the crossover scaling of p g using t x ~ z\. 
Indeed, numerical data collapse pretty well for z c = 5, 7, 
and 9, as shown in Fig. SJc). This can be explained by 
an intuitive argument. As discussed in the last part of 
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Sect. lIIIl the stochasticity inherent to toppling arrows re- 
veals only in the case when the unstable boundary site, 
say, (i,t) at the right boundary, has z c grains and its 
neighboring site at the next layer, say, (i + l,t + 1), is 
occupied by k grains. Assuming that the grain distribu- 
tion is uniform (to be confirmed), the probability of an 
occupied site having k grains is about z^ 1 . If the unsta- 
ble site has z c grains, it might have been occupied by k 
or k + 1 grains before receiving grains from the unstable 
neighboring site at the previous layer, whose probability 
is about z~ l . Then, the stochasticity in toppling arrows 
begins to appear when the expected number of such cases 
up to the layer t becomes finite, i.e., z~ 2 t ~ 0(1), which 
implies that t x ~ z 2 . 

In the aND at z c = 51, the value of a approaches 1 
in Fig. H](a), which may mislead to conclude no crossover 
for the larger z c values even if they are odd. We here 
pose the following question: How can we know whether or 
not any crossover behavior is observed at the sufficiently 
large layer? The value of a is not sufficient enough to 
answer the question. We need to investigate the grain 
distribution Q{z 1 t) in Eq. that contains the detailed 
information for metastable grain patterns. Due to the 
fact that the Q(0,t) = 1 — p g (t) becomes dominant as t 
increases for non-Abelian DSMs, we consider the rescaled 
grain distribution as 



Qr(z,t) = 



Q(z,t) 

i-Q(o,ty 



(21) 



where Q r (z, t) = for z G {1, 2, [z c /2] - 1}. 

We numerically find that Q r (z,t) is broad due to the 
stochastic nature for odd z c values. It is also found that 
it is peaked around at the value of [z c /2] due to the deter- 
ministic nature for even z c values. Figure HJd) shows the 
broadness of Q r (z,t) for z c = 51 and the peaked struc- 
ture of Q r (z, t) for z c = 50. This implies that the shape 
of the grain distribution is one of the good indicators to 
identify the universality class, as well as the grain density 
exponent in metastable patterns. 

The effect of the quenched randomness in the aND 
is compared to that of the annealed randomness in the 
pND. The scaling behaviors of p g (t) in the pND with 
p — i for various z c values exhibit qualitatively the same 
behaviors as those in the aND (not shown). Since the 
scaling behavior of p g (t) in the aND (a = i) at z c = 3 is 
quite different from that in the bND (a = 1) at the same 
z c value, one can ask if another crossover behavior exists 
between them. At each value of p close to 1 in the pND, 
we are able to observe the crossover behavior of the grain 
occupation density from a = 1 to ^ at some crossover 
layer, t x , as shown in Fig.[5](a). We numerically find that 
the value of t x scales as (1 — p) -1 / 2 , [see Fig.OJb)]. Such 
scaling is roughly explained by the following argument: 
the degree of stochasticity per toppling event is simply 
given by 1 — p. The total number of realizations up to 
the layer t is on the order of w(t)t ~ t 2 for t < t x . Thus, 
the stochasticity begins to appear when (1— p)t 2 ~ 0(1), 
which implies that t x ~ (1 — p)~~ x l 2 . 



Figure [UJa) shows the scaling behaviors of p e (t) in 
the bND for various z c values. The results imply the 
MF class with some logarithmic corrections to scaling 
for small z c values and without corrections for large z c 
values. For z c — 2 and 3, a cS (t) values are less than 1 
even in the scaling regime (not shown), but the possi- 
bility to approach 1 cannot be excluded. Such deviation 
affects other scaling behaviors. We also measure the ef- 
fective exponents of avalanche distributions and test the 
FSS collapse with MF values for z c = 2, 3, and 4 [see 
Fig. EJb)]. From our numerical observation, it is dif- 
ficult to determine whether the cases of z c — 2 and 3 
obey eventually the same FSS as the cases of z c > 4 or 
not. Figure |6jc) shows that Q r (z,t) is peaked around 
at [z c /2], consistent with the observation that avalanche 
boundaries behave ballistically, which implies again the 
MF class, independent of the parity of z c (> 4). 

Next, we investigate another kind of crossover behav- 
iors from the NS class to the MF class in the NS as 
z c increases up to 100. This is motivated by Liibeck's 
work [2l| that reported the crossover behavior between 
the undirected NS and the undirected cND for very large 
z c values (up to the order of 10 3 ). In the undirected NS, 
if to > z c grains topple at an unstable site, its neighbor- 
ing sites will receive grains of 0(to) on average with the 
standard deviation of 0(^/m). Thus, for sufficiently large 
z c values, the stochasticity in toppling rules is effectively 
suppressed and the system shows the deterministic uni- 
versality class. It was claimed that universality classes 
can be determined in the context of the shape of grain 
distributions: the grain distribution is broad for the undi- 
rected stochastic class, while it is narrow with multiple 
peaks for the undirected deterministic class. In the sim- 
ilar sense, one can expect essentially the same crossover 
from the NS class to the MF class in the directed case 
as z c increases to the sufficiently large value, i.e., by sup- 
pressing fluctuations in A,;±i it+ i. 

In Fig. E^a), we observe two kinds of crossover behav- 
iors of grain occupation densities in DSMs. One is from 
a = | (NS class) to 1 (MF class) as z c increases, as ex- 
pected. The other is from a = 1 to \ at a fixed z c value 
as t increases. Here the crossover layer, t Xl depends on 
the z c value. We systematically observe that the shape 
of Q r (z,t) at z c — 100 is much narrower than that at 
z c = 25 for any layer t [see Fig. [7[b)]. It implies the MF 
class for much larger z c values, which is hardly confirmed 
in numerical simulations. 

Finally, we discuss the crossover mechanism from the 
MF class to the NS class at a fixed z c value for the di- 
rected case. The mechanism is different from that for 
the undirected case. The numbers of toppled grains at 
unstable sites on the top layer t = are always z c . Then, 
Q(z,t — 1) has the normal distribution with the mean 
value of z c /2 and the standard deviation on the order of 
y/zZ, as well as an additional peak around z c — 1. The 
latter peak is due to the sites that have received grains 
both from the neighboring sites on the top layer, but 
are still stable. This peaked structure implies the MF 
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class. Inevitable fluctuations around two peaks in the 
Q(z, t) at small t propagate and accumulate through lay- 
ers, and finally result in the broad distribution without 
any distinguished peaks. If the value of z c was sufficiently 
large enough to exclude any fluctuations (approximately 
the cND), the MF behavior would be clearly observed 
without any crossover behaviors. In other words, for the 
finite z c values, there always exists crossovers from the 
MF class to the NS class at some finite crossover layer. 

V. SUMMARY 

We studied both non-Abelian stochastic (NS) and 
mean-field (MF) universality classes and crossover behav- 
iors between two classes in various non-Abelian directed 
sandpile models (DSMs) on a two-dimensional tilted lat- 
tice. The broken Abelian symmetry induces spatially 
long-range correlations in metastable patterns. The uni- 
versality class of avalanche dynamics can be identified by 
means of the grain density exponent in metastable pat- 
terns: a — I for the NS class and a = 1 for the MF 
class. Due to some ambiguity of the non-Abelian deter- 
ministic toppling rule against the given lattice structure, 
we considered the toppling bias for the last grain. We 
found that the parity of the threshold value, z c , might 
change the universality class. Using extensive numerical 
simulations, we confirmed that MF behaviors are well 
observed in the aND for even z c values and in the bND 
for any z c values as z c increases. For odd z c values, the 
toppling randomness in the aND and in the pND with 
p = \ becomes relevant, so that they belong to the NS 
class. 



In the aND and the pND with p = ^ for large odd z c 
values, we numerically observed crossover behaviors from 
a = 1 to J at some crossover layer, i x . The crossover 
layer is a function of the z c value and the partial bias, 
p, in the pND. Since the grain density exponent is not 
sufficient to determine the universality class under some 
circumstances, we suggest checking the broadness of the 
grain distribution in metastable patterns as well. Here 
the broadness originates from the stochastic nature in 
toppling rules. As in the undirected case, we also ob- 
served that the NS class undergoes the crossover to the 
MF class as z c increases. At a fixed z c value, another 
crossover is found from the MF class to the NS class at 
a finite value of t x . 

In conclusion, we emphasize that one should be very 
careful to set the toppling bias and the threshold value 
for some appropriate minimal model of self-organized 
criticality in order to explain experimental results. It 
is because the naive setting of parameters may mislead 
to either unexpected universality classes or crossover be- 
haviors. Finally, it would be worthwhile to extend our 
various tests towards some other lattice structures, to 
clarify the crossover scaling issues more analytically, and 
to test other types of non-Abelian DSMs, such as sticky 
sandpiles and sandpiles with bulk dissipation [261 ] . 
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